Dataset Match

We used the Treatment Outcomes Profile to join users of the C1 (n= 110,403; users=85,603) that were evaluated with the Treatment Outcome Profile (TOP) in the same date of Admission to entries in C1 datasets.


###casos duplicados en fecha de aplicación
#1     1 97779
#2     2  2483
#3     3   149
#4     4    21
#5     5     3
#no hay casos que tengan más de una fecha de admisión Y UNA BD MÁS RECIENTE
#CONS_TOP_df_dup_ENE_2020_prev9%>% dplyr::mutate(concat=paste0(hash_key,"_",as.character(fech_ap_top)))%>% group_by(concat) %>% dplyr::mutate(n_hash_fech_ap=n(),n_dis_ano=n_distinct(ano_bd))%>% dplyr::filter(n_hash_fech_ap>1,n_dis_ano>1)

invisible(
CONS_TOP_df_dup_ENE_2020_prev9%>% dplyr::mutate(concat=paste0(hash_key,"_",as.character(fech_ap_top)))%>% group_by(concat) %>% dplyr::mutate(n_hash_fech_ap=n(),n_dis_ano=n_distinct(ano_bd))%>% dplyr::filter(n_hash_fech_ap>1,n_dis_ano>1)
)

CONS_TOP_df_dup_ENE_2020_prev9_match<-
CONS_TOP_df_dup_ENE_2020_prev9%>% dplyr::mutate(concat=paste0(hash_key,"_",as.character(fech_ap_top)))%>% group_by(concat) %>% dplyr::mutate(n_hash_fech_ap=n(),n_dis_ano=n_distinct(ano_bd))%>% dplyr::filter(n_hash_fech_ap==1)%>%
  dplyr::select(hash_key,concat,total_oh,dosis_oh,total_thc,dosis_thc,total_pbc,dosis_pbc,total_coc,dosis_coc,total_bzd,dosis_bzd,total_otra,dosis_otra,hurto,robo,venta_drogas,rina,total_vif,otro,total_transgresion,salud_psicologica,total_trabajo,total_educacion,salud_fisica,lugar_vivir,vivienda,calidad_vida,etapa_del_tratamiento)

CONS_C1_df_dup_JUL_2020_match_top<-
  CONS_C1_df_dup_JUL_2020%>%
  dplyr::mutate(concat=paste0(hash_key,"_",as.character(fech_ing)))%>%
  dplyr::left_join(CONS_TOP_df_dup_ENE_2020_prev9_match,by="concat", suffix=c("","_TOP"))%>%
  dplyr::filter(!is.na(compromiso_biopsicosocial), etapa_del_tratamiento=="Inicio Tratamiento")%>%
  dplyr::select(-hash_key_TOP)%>%
  dplyr::filter_at(vars(total_oh,dosis_oh,total_thc,dosis_thc,total_pbc,dosis_pbc,total_coc,dosis_coc,total_bzd,dosis_bzd,total_otra,dosis_otra,hurto,robo,venta_drogas,rina,total_vif,otro,total_transgresion,salud_psicologica,total_trabajo,total_educacion,salud_fisica,lugar_vivir,vivienda,calidad_vida,etapa_del_tratamiento), all_vars(!is.na(.)))

    users_top_only_one_app_match<-CONS_TOP_df_dup_ENE_2020_prev9%>% dplyr::mutate(concat=paste0(hash_key,"_",as.character(fech_ap_top)))%>% group_by(concat) %>% dplyr::mutate(n_hash_fech_ap=n(),n_dis_ano=n_distinct(ano_bd))%>% dplyr::filter(n_hash_fech_ap==1)%>% dplyr::ungroup()%>% dplyr::distinct(hash_key)%>% nrow()


#c("total_oh","dosis_oh","total_thc","dosis_thc","total_pbc","dosis_pbc","total_coc","dosis_coc","total_bzd","dosis_bzd","total_otra","dosis_otra","hurto","robo","venta_drogas","rina","total_vif","otro","total_transgresion","salud_psicologica","total_trabajo","total_educacion","salud_fisica","lugar_vivir","vivienda","calidad_vida")

# etapa_del_tratamiento      n      percent valid_percent
#                Egreso      4 3.470114e-05  0.0005719188
#    Inicio Tratamiento   6948 6.027587e-02  0.9934229339
#  Seguimiento 12 meses      0 0.000000e+00  0.0000000000
#  Seguimiento 15 meses      0 0.000000e+00  0.0000000000
#   Seguimiento 3 meses     34 2.949597e-04  0.0048613097
#   Seguimiento 6 meses      8 6.940227e-05  0.0011438376
#   Seguimiento 9 meses      0 0.000000e+00  0.0000000000
#                  <NA> 108276 9.393251e-01            NA


Of the total entries in TOP dataset (n= 110,403), we did not considered 2,845 entries in TOP dataset had more than one application in the same date (3%). We selected 6,390 cases that had no differences between the date of application and the date of admission to treatment. Additionally, these applications must had been administered at admission (excluding the 0,6% of applications that were administered at other stage of treatment). Finally, we discarded cases with null values in the variables of interest (230 entries).


Explore variables

First, values of doses over 50 were discarded and declared missing, due to the great amount of outliers.


We defined one outcome variable to indicate early discharge (abandono_temp) that should vary depending on biopsychosocial compromise. We also collapsed variables related to transgression to norms (transg_norm) into a continuous variable, by counting actions committed in the last 4 weeks. Additionally, we added the variable salud_mean to get the mean of the scores in physical health, quality of life and psychological health. To capture the substance use behaviors, we defined the average pattern of drug use in alcohol, marijuana, cocaine paste base, cocaine sedatives and tranquilizers, and other substances, in terms of scores of quantity in the last four weeks (total_mean). Equally, we defined a score of the average dose of the mean reported by users on each substance. At last, we created the mean of scores of days of attendance to jobs or educational institutions in the last four weeks (trabajo_edu_mean).


Furthermore, we were interested in co-occurrent physical and psychiatric conditions among substance users. Consequently, we generated a variable called dg_cie_10, composed of users that at least had been diagnosed with a psychiatric condition or professionals had been studying a possible diagnosis. Also, we included dg_trs_fisico, which categorizes between users without a diagnosed physical condition (0), in study (1), and those diagnosed with a physical condition (2).


To capture profiles of consumption, we generated the following variables: oh identifies cases that used alcohol at least one day in the last four weeks. This same criteria was applied to marijuana (thc), cocaine paste base (pbc), cocaine (coc), benzodiazepines (bzd), and other drugs (otra). Users that reported at least two substances according to the aforementioned rule were categorized as polydrug users (polycons).


Lastly, the mean of the scores in many variables were dichotomized between users that had less than the 50% of users, compared to people that showed more scores than the 50% of the population. These were defined as Higher Frequency of Days using substances in the last four weeks (mayor_total), Higher Average Reported Daily Dose of Substances (mayor_dosis), Lower Health-Related Self-Reports (menor_salud), and Lower Days of Attendance to Work and Educational Centers in the last four weeks (menor_trabajo_edu). In a similar manner, transg_norm_bin identify users that reported at least one behavior that transgress the norm.


Warning in sorted_count(x): Variable contains value(s) of "" that have been
converted to "empty".
Table 1. Summary of Variables
Type Variable Missing Complete Rate Ordered No. Unique Top Counts Mean (logical) Counts (logical) Mean SD Min Perc. 25 Median Perc. 75 Max Histogram
factor compromiso_biopsicosocial 0 1.0000000 FALSE 3 2-M: 3593, 3-S: 1976, 1-L: 810 NA NA NA NA NA NA NA NA NA NA
factor dg_cie_10 0 1.0000000 FALSE 2 1: 3884, 0: 2495 NA NA NA NA NA NA NA NA NA NA
factor dg_trs_fisico 0 1.0000000 FALSE 3 1: 3648, 0: 2065, 2: 666 NA NA NA NA NA NA NA NA NA NA
factor abandono_temp 0 1.0000000 FALSE 2 0: 5163, 1: 1216 NA NA NA NA NA NA NA NA NA NA
factor origen_ingreso_mod 0 1.0000000 FALSE 5 Con: 3441, Sec: 1577, Der: 557, Sec: 516 NA NA NA NA NA NA NA NA NA NA
factor motivodeegreso_mod_imp 0 1.0000000 FALSE 7 Aba: 2047, Alt: 1567, Aba: 1216, Der: 645 NA NA NA NA NA NA NA NA NA NA
factor dg_trs_cons_sus_or 0 1.0000000 FALSE 2 Dep: 4454, Con: 1925 NA NA NA NA NA NA NA NA NA NA
factor diagnostico_trs_fisico 0 1.0000000 FALSE 25 En : 3648, Sin: 2065, Otr: 168, Tra: 115 NA NA NA NA NA NA NA NA NA NA
factor otros_probl_at_sm_or 1,437 0.7747296 FALSE 14 Sin: 2564, Vio: 1587, Otr: 625, Abu: 96 NA NA NA NA NA NA NA NA NA NA
factor sus_principal_mod 0 1.0000000 FALSE 5 Alc: 2468, Pas: 2204, Coc: 1305, Mar: 323 NA NA NA NA NA NA NA NA NA NA
factor otras_sus1_mod 2,102 0.6704813 FALSE 5 Alc: 1752, Mar: 1266, Coc: 732, Pas: 350 NA NA NA NA NA NA NA NA NA NA
factor otras_sus2_mod 4,109 0.3558551 FALSE 5 Mar: 861, Alc: 657, Coc: 434, Pas: 164 NA NA NA NA NA NA NA NA NA NA
factor otras_sus3_mod 4,109 0.3558551 FALSE 5 Mar: 861, Alc: 657, Coc: 434, Pas: 164 NA NA NA NA NA NA NA NA NA NA
factor alta_adm 0 1.0000000 FALSE 2 0: 5828, 1: 551 NA NA NA NA NA NA NA NA NA NA
logical hurto 0 1.0000000 NA NA NA 0.1067565 FAL: 5698, TRU: 681 NA NA NA NA NA NA NA NA
logical robo 0 1.0000000 NA NA NA 0.0608246 FAL: 5991, TRU: 388 NA NA NA NA NA NA NA NA
logical venta_drogas 0 1.0000000 NA NA NA 0.0443643 FAL: 6096, TRU: 283 NA NA NA NA NA NA NA NA
logical rina 0 1.0000000 NA NA NA 0.1250980 FAL: 5581, TRU: 798 NA NA NA NA NA NA NA NA
logical otro 0 1.0000000 NA NA NA 0.0335476 FAL: 6165, TRU: 214 NA NA NA NA NA NA NA NA
logical lugar_vivir 0 1.0000000 NA NA NA 0.9211475 TRU: 5876, FAL: 503 NA NA NA NA NA NA NA NA
logical vivienda 0 1.0000000 NA NA NA 0.9230287 TRU: 5888, FAL: 491 NA NA NA NA NA NA NA NA
logical vif 0 1.0000000 NA NA NA 0.2022261 FAL: 5089, TRU: 1290 NA NA NA NA NA NA NA NA
logical transg_norm_bin 0 1.0000000 NA NA NA 0.3563254 FAL: 4106, TRU: 2273 NA NA NA NA NA NA NA NA
logical menor_vivienda 0 1.0000000 NA NA NA 0.0987616 FAL: 5749, TRU: 630 NA NA NA NA NA NA NA NA
logical menor_salud 0 1.0000000 NA NA NA 0.4400376 FAL: 3572, TRU: 2807 NA NA NA NA NA NA NA NA
logical menor_trabajo_edu 0 1.0000000 NA NA NA 0.4782881 FAL: 3328, TRU: 3051 NA NA NA NA NA NA NA NA
logical mayor_total 0 1.0000000 NA NA NA 0.5063490 TRU: 3230, FAL: 3149 NA NA NA NA NA NA NA NA
logical mayor_dosis 0 1.0000000 NA NA NA 0.5099545 TRU: 3253, FAL: 3126 NA NA NA NA NA NA NA NA
logical polycons 0 1.0000000 NA NA NA 0.8935570 TRU: 5700, FAL: 679 NA NA NA NA NA NA NA NA
logical mayor_comp_biopsicsoc 0 1.0000000 NA NA NA 0.3097664 FAL: 4403, TRU: 1976 NA NA NA NA NA NA NA NA
logical menor_comp_biopsicsoc 0 1.0000000 NA NA NA 0.1269792 FAL: 5569, TRU: 810 NA NA NA NA NA NA NA NA
numeric total_oh 0 1.0000000 NA NA NA NA NA 7.6126352 8.7881215 0.0000000 0.0000000 4.0000000 12.0000000 28.0000000 <U+2587><U+2582><U+2582><U+2581><U+2582>
numeric dosis_oh 0 1.0000000 NA NA NA NA NA 0.8424518 2.1179793 0.0000000 0.0000000 0.0000000 1.0000000 39.0000000 <U+2587><U+2581><U+2581><U+2581><U+2581>
numeric total_thc 0 1.0000000 NA NA NA NA NA 4.6778492 8.9716897 0.0000000 0.0000000 0.0000000 4.0000000 28.0000000 <U+2587><U+2581><U+2581><U+2581><U+2581>
numeric dosis_thc 0 1.0000000 NA NA NA NA NA 0.8424518 2.1179793 0.0000000 0.0000000 0.0000000 1.0000000 39.0000000 <U+2587><U+2581><U+2581><U+2581><U+2581>
numeric total_pbc 0 1.0000000 NA NA NA NA NA 5.0953127 9.1610794 0.0000000 0.0000000 0.0000000 6.0000000 28.0000000 <U+2587><U+2581><U+2581><U+2581><U+2581>
numeric dosis_pbc 0 1.0000000 NA NA NA NA NA 1.1705596 2.9942045 0.0000000 0.0000000 0.0000000 1.0000000 40.0000000 <U+2587><U+2581><U+2581><U+2581><U+2581>
numeric total_coc 0 1.0000000 NA NA NA NA NA 2.6963474 6.1555415 0.0000000 0.0000000 0.0000000 1.0000000 28.0000000 <U+2587><U+2581><U+2581><U+2581><U+2581>
numeric dosis_coc 0 1.0000000 NA NA NA NA NA 0.2688509 1.7629119 0.0000000 0.0000000 0.0000000 0.0000000 40.0000000 <U+2587><U+2581><U+2581><U+2581><U+2581>
numeric total_bzd 0 1.0000000 NA NA NA NA NA 0.8076501 4.0691149 0.0000000 0.0000000 0.0000000 0.0000000 28.0000000 <U+2587><U+2581><U+2581><U+2581><U+2581>
numeric dosis_bzd 0 1.0000000 NA NA NA NA NA 0.4514814 2.8107924 0.0000000 0.0000000 0.0000000 0.0000000 40.0000000 <U+2587><U+2581><U+2581><U+2581><U+2581>
numeric total_otra 0 1.0000000 NA NA NA NA NA 0.9904374 4.9582202 0.0000000 0.0000000 0.0000000 0.0000000 28.0000000 <U+2587><U+2581><U+2581><U+2581><U+2581>
numeric dosis_otra 0 1.0000000 NA NA NA NA NA 0.4514814 2.8107924 0.0000000 0.0000000 0.0000000 0.0000000 40.0000000 <U+2587><U+2581><U+2581><U+2581><U+2581>
numeric total_vif 0 1.0000000 NA NA NA NA NA 1.3500549 4.5085774 0.0000000 0.0000000 0.0000000 0.0000000 28.0000000 <U+2587><U+2581><U+2581><U+2581><U+2581>
numeric total_transgresion 0 1.0000000 NA NA NA NA NA 0.3705910 0.7909335 0.0000000 0.0000000 0.0000000 0.0000000 5.0000000 <U+2587><U+2581><U+2581><U+2581><U+2581>
numeric salud_psicologica 0 1.0000000 NA NA NA NA NA 9.2128860 4.6302091 0.0000000 5.0000000 10.0000000 12.0000000 20.0000000 <U+2583><U+2586><U+2587><U+2585><U+2581>
numeric total_trabajo 0 1.0000000 NA NA NA NA NA 10.9487380 10.5422163 0.0000000 0.0000000 10.0000000 20.0000000 28.0000000 <U+2587><U+2581><U+2582><U+2583><U+2583>
numeric total_educacion 0 1.0000000 NA NA NA NA NA 0.3099232 2.2582653 0.0000000 0.0000000 0.0000000 0.0000000 28.0000000 <U+2587><U+2581><U+2581><U+2581><U+2581>
numeric salud_fisica 0 1.0000000 NA NA NA NA NA 11.1327794 4.9691828 0.0000000 8.0000000 10.0000000 15.0000000 20.0000000 <U+2582><U+2585><U+2587><U+2586><U+2583>
numeric calidad_vida 0 1.0000000 NA NA NA NA NA 10.4111930 5.1114644 0.0000000 7.0000000 10.0000000 15.0000000 20.0000000 <U+2583><U+2585><U+2587><U+2585><U+2583>
numeric sin_otros_prob_sm 0 1.0000000 NA NA NA NA NA 0.4021006 0.4903605 0.0000000 0.0000000 0.0000000 1.0000000 1.0000000 <U+2587><U+2581><U+2581><U+2581><U+2586>
numeric dg_cons_dep 0 1.0000000 NA NA NA NA NA 0.6982286 0.4590626 0.0000000 0.0000000 1.0000000 1.0000000 1.0000000 <U+2583><U+2581><U+2581><U+2581><U+2587>
numeric dg_trs_psiq_cie_10_or 0 1.0000000 NA NA NA NA NA 1.0039191 0.8970966 0.0000000 0.0000000 1.0000000 2.0000000 2.0000000 <U+2587><U+2581><U+2583><U+2581><U+2587>
numeric x2_dg_trs_psiq_cie_10_or 0 1.0000000 NA NA NA NA NA 0.0852798 0.3780609 0.0000000 0.0000000 0.0000000 0.0000000 2.0000000 <U+2587><U+2581><U+2581><U+2581><U+2581>
numeric x3_dg_trs_psiq_cie_10_or 0 1.0000000 NA NA NA NA NA 0.0227308 0.1922389 0.0000000 0.0000000 0.0000000 0.0000000 2.0000000 <U+2587><U+2581><U+2581><U+2581><U+2581>
numeric x4_dg_trs_psiq_cie_10_or 0 1.0000000 NA NA NA NA NA 0.0006271 0.0354107 0.0000000 0.0000000 0.0000000 0.0000000 2.0000000 <U+2587><U+2581><U+2581><U+2581><U+2581>
numeric x5_dg_trs_psiq_cie_10_or 0 1.0000000 NA NA NA NA NA 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 <U+2581><U+2581><U+2587><U+2581><U+2581>
numeric x6_dg_trs_psiq_cie_10_or 0 1.0000000 NA NA NA NA NA 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 <U+2581><U+2581><U+2587><U+2581><U+2581>
numeric transg_norm 0 1.0000000 NA NA NA NA NA 0.5728171 0.9533594 0.0000000 0.0000000 0.0000000 1.0000000 6.0000000 <U+2587><U+2581><U+2581><U+2581><U+2581>
numeric salud_mean 0 1.0000000 NA NA NA NA NA 10.2522861 4.1175032 0.0000000 7.3333333 10.0000000 13.3333333 20.0000000 <U+2582><U+2585><U+2587><U+2586><U+2582>
numeric dosis_mean 0 1.0000000 NA NA NA NA NA 0.6712128 1.3873650 0.0000000 0.0000000 0.1666667 0.6666667 19.0000000 <U+2587><U+2581><U+2581><U+2581><U+2581>
numeric total_mean 0 1.0000000 NA NA NA NA NA 3.6467053 3.3507752 0.0000000 0.8333333 2.8333333 5.3333333 23.3333333 <U+2587><U+2583><U+2581><U+2581><U+2581>
numeric vivienda_mean 0 1.0000000 NA NA NA NA NA 4.6104405 1.2392178 0.0000000 5.0000000 5.0000000 5.0000000 5.0000000 <U+2581><U+2581><U+2581><U+2581><U+2587>
numeric trabajo_edu_mean 0 1.0000000 NA NA NA NA NA 5.6293306 5.3284457 0.0000000 0.0000000 5.0000000 10.0000000 24.0000000 <U+2587><U+2582><U+2586><U+2581><U+2581>
numeric total_cie_10 0 1.0000000 NA NA NA NA NA 1.1125568 1.0945562 0.0000000 0.0000000 1.0000000 2.0000000 8.0000000 <U+2587><U+2585><U+2581><U+2581><U+2581>
numeric salud_mean_z 0 1.0000000 NA NA NA NA NA 0.0013466 0.9996963 -2.4878250 -0.7073514 -0.0599064 0.7493998 2.3680123 <U+2582><U+2585><U+2587><U+2586><U+2582>
numeric dosis_mean_z 0 1.0000000 NA NA NA NA NA 0.0000000 1.0000000 -0.4838041 -0.4838041 -0.3636723 -0.0032768 13.2112223 <U+2587><U+2581><U+2581><U+2581><U+2581>
numeric total_mean_z 0 1.0000000 NA NA NA NA NA -0.0007103 0.9992102 -1.0881676 -0.8396653 -0.2432600 0.5022467 5.8698950 <U+2587><U+2583><U+2581><U+2581><U+2581>
numeric vivienda_mean_z 0 1.0000000 NA NA NA NA NA -0.0002256 1.0005570 -3.7227419 0.3143087 0.3143087 0.3143087 0.3143087 <U+2581><U+2581><U+2581><U+2581><U+2587>
numeric total_vif_z 0 1.0000000 NA NA NA NA NA -0.0009058 0.9979097 -0.2997213 -0.2997213 -0.2997213 -0.2997213 5.8976821 <U+2587><U+2581><U+2581><U+2581><U+2581>
numeric oh 0 1.0000000 NA NA NA NA NA 0.7239379 0.4470829 0.0000000 0.0000000 1.0000000 1.0000000 1.0000000 <U+2583><U+2581><U+2581><U+2581><U+2587>
numeric thc 0 1.0000000 NA NA NA NA NA 0.3411193 0.4741225 0.0000000 0.0000000 0.0000000 1.0000000 1.0000000 <U+2587><U+2581><U+2581><U+2581><U+2585>
numeric pbc 0 1.0000000 NA NA NA NA NA 0.3461358 0.4757744 0.0000000 0.0000000 0.0000000 1.0000000 1.0000000 <U+2587><U+2581><U+2581><U+2581><U+2585>
numeric coc 0 1.0000000 NA NA NA NA NA 0.3049067 0.4604040 0.0000000 0.0000000 0.0000000 1.0000000 1.0000000 <U+2587><U+2581><U+2581><U+2581><U+2583>
numeric bzd 0 1.0000000 NA NA NA NA NA 0.0647437 0.2460923 0.0000000 0.0000000 0.0000000 0.0000000 1.0000000 <U+2587><U+2581><U+2581><U+2581><U+2581>
numeric otra 0 1.0000000 NA NA NA NA NA 0.0465590 0.2107089 0.0000000 0.0000000 0.0000000 0.0000000 1.0000000 <U+2587><U+2581><U+2581><U+2581><U+2581>
#knitr::include_graphics("G:/Mi unidad/Alvacast/SISTRAT 2019 (github)/SUD_CL/Figures/Figure_Duplicates.svg")
#foreign::write.dta(CONS_C1_df_dup_JUL_2020_match_top_sel, "G:/Mi unidad/Alvacast/SISTRAT 2019 (github)/c1_top.dta")

#CONS_C1_df_dup_JUL_2020 CONS_C1_df_dup_JUL_2020_match_top CONS_TOP_df_dup_ENE_2020_prev9_match CONS_TOP_df_dup_ENE_2020_prev9
#sólo match
comb_datasets_a_n<-CONS_C1_df_dup_JUL_2020%>%
    dplyr::mutate(concat=paste0(hash_key,"_",as.character(fech_ing)))%>%
    dplyr::left_join(CONS_TOP_df_dup_ENE_2020_prev9_match,by="concat", suffix=c("","_TOP"))%>%
    dplyr::filter(!is.na(hash_key_TOP))%>%
    dplyr::select(-hash_key_TOP)%>% 
    nrow()
comb_datasets_a_users<-CONS_C1_df_dup_JUL_2020%>%
    dplyr::mutate(concat=paste0(hash_key,"_",as.character(fech_ing)))%>%
    dplyr::left_join(CONS_TOP_df_dup_ENE_2020_prev9_match,by="concat", suffix=c("","_TOP"))%>%
    dplyr::filter(!is.na(hash_key_TOP))%>%
    dplyr::select(-hash_key_TOP)%>% 
    dplyr::distinct(hash_key)%>% nrow()
#sólo los que tienen compromiso biopsicosocial
comb_datasets_b_n<-CONS_C1_df_dup_JUL_2020%>%
    dplyr::mutate(concat=paste0(hash_key,"_",as.character(fech_ing)))%>%
    dplyr::left_join(CONS_TOP_df_dup_ENE_2020_prev9_match,by="concat", suffix=c("","_TOP"))%>%
    dplyr::filter(!is.na(compromiso_biopsicosocial), !is.na(hash_key_TOP))%>%
    dplyr::select(-hash_key_TOP)%>% 
    nrow()
comb_datasets_b_users<-CONS_C1_df_dup_JUL_2020%>%
    dplyr::mutate(concat=paste0(hash_key,"_",as.character(fech_ing)))%>%
    dplyr::left_join(CONS_TOP_df_dup_ENE_2020_prev9_match,by="concat", suffix=c("","_TOP"))%>%
    dplyr::filter(!is.na(compromiso_biopsicosocial), !is.na(hash_key_TOP))%>%
    dplyr::select(-hash_key_TOP)%>% 
    dplyr::distinct(hash_key)%>% nrow()
#sólo los que tienen etapa de tratamiento de inicio de tratamiento
comb_datasets_c_n<-CONS_C1_df_dup_JUL_2020%>%
    dplyr::mutate(concat=paste0(hash_key,"_",as.character(fech_ing)))%>%
    dplyr::left_join(CONS_TOP_df_dup_ENE_2020_prev9_match,by="concat", suffix=c("","_TOP"))%>%
    dplyr::filter(!is.na(compromiso_biopsicosocial), !is.na(hash_key_TOP), etapa_del_tratamiento=="Inicio Tratamiento")%>%
    dplyr::select(-hash_key_TOP)%>% 
    nrow()
comb_datasets_c_users<-CONS_C1_df_dup_JUL_2020%>%
    dplyr::mutate(concat=paste0(hash_key,"_",as.character(fech_ing)))%>%
    dplyr::left_join(CONS_TOP_df_dup_ENE_2020_prev9_match,by="concat", suffix=c("","_TOP"))%>%
    dplyr::filter(!is.na(compromiso_biopsicosocial), !is.na(hash_key_TOP), etapa_del_tratamiento=="Inicio Tratamiento")%>%
    dplyr::select(-hash_key_TOP)%>% 
    dplyr::distinct(hash_key)%>% nrow()

#https://stackoverflow.com/questions/46750364/diagrammer-and-graphviz
#https://mikeyharper.uk/flowcharts-in-r-using-diagrammer/
#http://blog.nguyenvq.com/blog/2012/05/29/better-decision-tree-graphics-for-rpart-via-party-and-partykit/
#http://blog.nguyenvq.com/blog/2014/01/17/skeleton-to-create-fast-automatic-tree-diagrams-using-r-and-graphviz/
#https://cran.r-project.org/web/packages/DiagrammeR/vignettes/graphviz-mermaid.html
#https://stackoverflow.com/questions/39133058/how-to-use-graphviz-graphs-in-diagrammer-for-r
#https://subscription.packtpub.com/book/big_data_and_business_intelligence/9781789802566/1/ch01lvl1sec21/creating-diagrams-via-the-diagrammer-package
#https://justlegal.be/2019/05/using-flowcharts-to-display-legal-procedures/
#
#
library(DiagrammeR) #⋉
grViz("digraph flowchart {
      # node definitions with substituted label text
      node [fontname = Helvetica, shape = rectangle,fontsize = 10]        
      tab1 [label = '@@1']
      tab2 [label = '@@2']
      tab3 [label = '@@3']
      tab4 [label = '@@4']
      tab5 [label = '@@5']
      tab6 [label = '@@6']
      tab7 [label = '@@7']
      blank [label = '', width = 0.001, height = 0.001]

      # edge definitions with the node IDs
      tab1 -> tab2;
      blank -> tab3[label='',dir = none, shape=none, color = 'white',fontcolor = black, width=0, height=0];
      {rank=same; tab2 -> tab3 [label='Left join',fontsize = 11, dir=none]}; #⋉
      tab3 -> tab4 [label='  Result of the left join of the dataset',fontsize = 9];
      tab4 -> tab5 [label='  Only rows with available data on Biopsychosocial Compromise',fontsize = 9];
      tab5 -> tab6 [label='  Only strict matches with applications at the stage of admission',fontsize = 9];
      tab6 -> tab7 [label='  Only rows with available data on TOP scores and Diagnostic of CIE-10',fontsize = 9];
      }
Only rows with available data on TOP scores

      [1]:  paste0('TOP Dataset \\n(n = ', formatC(nrow(CONS_TOP_df_dup_ENE_2020_prev9), format='f', big.mark=',', digits=0), ';\\n users: ',formatC(CONS_TOP_df_dup_ENE_2020_prev9%>% dplyr::distinct(hash_key)%>% nrow(), format='f', big.mark=',', digits=0),')')
      [2]:  paste0('Only applications w/ only one\\n application in the same date \\n(n = ', formatC(nrow(CONS_TOP_df_dup_ENE_2020_prev9_match), format='f', big.mark=',', digits=0), ';\\n users:',formatC(users_top_only_one_app_match, format='f', big.mark=',', digits=0),')')
      [3]:  paste0('Only applications w/ only one\\n application in the same date \\n(n = ', formatC(nrow(CONS_C1_df_dup_JUL_2020), format='f', big.mark=',', digits=0), ';\\n users:',formatC(CONS_C1_df_dup_JUL_2020%>% dplyr::distinct(hash_key)%>% nrow(), format='f', big.mark=',', digits=0),')')
      [4]:  paste0('Dataset \\n(n = ',formatC(comb_datasets_a_n,format='f', big.mark=',', digits=0),'\\nusers =',formatC(comb_datasets_a_users,format='f', big.mark=',', digits=0),')')
      [5]:  paste0('Dataset \\n(n = ',formatC(comb_datasets_b_n,format='f', big.mark=',', digits=0),'\\nusers =',formatC(comb_datasets_b_users,format='f', big.mark=',', digits=0),')')
      [6]:  paste0('Dataset \\n(n = ',formatC(comb_datasets_c_n,format='f', big.mark=',', digits=0),'\\nusers =',formatC(comb_datasets_c_users,format='f', big.mark=',', digits=0),')')
      [7]:  paste0('Final Sample \\n(n = ', formatC(nrow(CONS_C1_df_dup_JUL_2020_match_top_sel), format='f', big.mark=',', digits=0), ';\\n users: ',formatC(CONS_C1_df_dup_JUL_2020_match_top_sel%>% dplyr::distinct(hash_key)%>% nrow(), format='f', big.mark=',', digits=0),')')
      ")


To get an idea of the distribution of the biopsychosocial compromise among the selected sample, we generated a pie chart.


library(plotly)

CONS_C1_df_dup_JUL_2020_match_top_sel %>%
     dplyr::group_by(compromiso_biopsicosocial)%>%
      dplyr::mutate(compromiso_biopsicosocial=case_when(compromiso_biopsicosocial=="2-Moderado"~"Moderate",
                                                        compromiso_biopsicosocial=="1-Leve"~"Mild",
                                                        TRUE~"Severe"))%>%
       dplyr::summarise(n=n(),percent=round(n/sum(n),2))%>% 
     dplyr::ungroup()%>%
plot_ly( labels = ~compromiso_biopsicosocial, values = ~n, type = 'pie',hole = 0.6,textposition = 'outside',textinfo = 'label+percent', marker = list(colors = c("yellowgreen","orange","tomato"))) %>%
  layout(xaxis = list(showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE),
         yaxis = list(showgrid = FALSE, zeroline = FALSE, showticklabels = FALSE))

Figure 2. Biopsychosocial Compromise

We generated a correlation plot which wraps different type of variables (continuous, polytomous, dichotomous)1. We selected mostly the original variables from the instrument.


# It first identifies the types of variables, organizes them by type (continuous, polytomous, dichotomous), calls the appropriate correlation function, and then binds the resulting matrices together.D
mix_cor<- psych::mixedCor(data=CONS_C1_df_dup_JUL_2020_match_top_sel%>%dplyr::mutate(compromiso_biopsicosocial=as.numeric(compromiso_biopsicosocial)),c=2:20, #total_oh,dosis_oh total_thc dosis_thc total_pbc dosis_pbc total_coc dosis_coc total_bzd dosis_bzd  total_otra dosis_otra total_vif  total_transgresion salud_psicologica total_trabajo total_educacion salud_fisica calidad_vida 
         p=1,#compromiso_biopsicosocial,
         d=21:27,#, hurto robo  venta_drogas rina  otro lugar_vivir vivienda
         smooth=TRUE,
         correct=.5,
         global=TRUE,
         ncat=8,
         use="pairwise",
         method="kendall",
         weight=NULL)
# Getting the correlation matrix of the dataset.
comp_biopsicosoc_mat <- mix_cor$rho

# Explore the correlation structure of the dataset.
library(ggcorrplot)
ggcorrplot(comp_biopsicosoc_mat,
          ggtheme = ggplot2::theme_gray,
          colors = c("#E46726", "white", "#6D9EC1"),
           tl.cex=8)
Figure 3.Matrix Of Correlations Between Original Variables

Figure 3.Matrix Of Correlations Between Original Variables

#subscript out of bounds

#mixedCor(data=CONS_C1_df_dup_JUL_2020_match_top_sel,c=c("total_oh","dosis_oh", "total_thc", "dosis_thc", "total_pbc", "dosis_pbc", "total_coc", "dosis_coc", "total_bzd", "dosis_bzd", "total_otra", "dosis_otra", "total_vif",  "total_transgresion", "salud_psicologica", "total_trabajo", "total_educacion", "salud_fisica", "calidad_vida"), p="compromiso_biopsicosocial",d=c("hurto", "robo",  "venta_drogas", "rina",  "otro", "lugar_vivir", "vivienda"), smooth=TRUE,correct=.5,global=TRUE,ncat=8,use="pairwise",method="pearson",weight=NULL)

#subscript out of bounds
#set.seed(123)
#par(mar=c(1,1,1,1))
#comp_biopsicosoc_mat_ci<-mycor.ci(x=CONS_C1_df_dup_JUL_2020_match_top_sel%>%dplyr::mutate(compromiso_biopsicosocial=as.numeric(compromiso_biopsicoso#cial)), 
#                                  method="spearman", 
#                                  poly=TRUE, 
#                                  cvars=2:20, 
#                                  pvars=1,
#                                  dvars=21:27, 
#                                   n.iter=1000,
#                                  plot = F)
#
#cor_pmat(): para ver la sig estadística de las correlaciones


Also we included an heterogeneus matrix of correlations, consisting of polychoric, biserial and pearson correlations.2 We mostly selected transformed variables in the following matrix.


require(polycor)
#Computes a heterogenous correlation matrix, consisting of Pearson product-moment correlations between numeric variables, polyserial correlations between numeric and ordinal variables, and polychoric correlations between ordinal variables.
hetcor_mat<-hetcor(CONS_C1_df_dup_JUL_2020_match_top_sel[c("compromiso_biopsicosocial","lugar_vivir","vivienda","transg_norm","salud_mean","dosis_mean","total_mean","trabajo_edu_mean","vivienda_mean","polycons", "menor_salud", "menor_trabajo_edu", "mayor_dosis", "mayor_total","transg_norm_bin","dg_cons_dep", "sin_otros_prob_sm","abandono_temp","alta_adm","dg_cie_10")], ML = T, std.err = T, use="complete.obs", bins=4, pd=TRUE)
#CONS_C1_df_dup_JUL_2020_match_top_sel%>% select(contains("dosis"))%>% hist()

ggcorrplot(hetcor_mat$correlations,
          ggtheme = ggplot2::theme_gray,
          colors = c("#E46726", "white", "#6D9EC1"), tl.cex=8)
Figure 4.Matrix Of Correlations Between Transformed Variables

Figure 4.Matrix Of Correlations Between Transformed Variables

#pcor(c("a", "b", "x", "y", "z"), var(mydata))
# partial corr between a and b controlling for x, y, z


Variables of TOP depending on Biopsychosocial Compromise


library(reshape2)  # for melt() function

Attaching package: 'reshape2'
The following objects are masked from 'package:data.table':

    dcast, melt
The following object is masked from 'package:tidyr':

    smiths
library(ggplot2)
# First we need to restructure the data into long format:
ocdMelt <- melt(CONS_C1_df_dup_JUL_2020_match_top_sel[c('compromiso_biopsicosocial','hurto',"robo","venta_drogas","rina","otro","lugar_vivir","vivienda","polycons","vif","transg_norm_bin")], id=c('compromiso_biopsicosocial'))
names(ocdMelt) <- c('Group', 'Outcome_Measure', 'Frequency')
ocdMelt <-ocdMelt%>% dplyr::mutate(Frequency=factor(Frequency))%>% dplyr::group_by(Outcome_Measure,Group,Frequency)%>% dplyr::arrange(Outcome_Measure, Group,Frequency)%>% dplyr::add_tally()%>% distinct(.,.keep_all=T)%>% dplyr::ungroup()%>% dplyr::group_by(Group, Outcome_Measure)%>% dplyr::mutate(perc=n/sum(n))%>% dplyr::ungroup()

# plot
p5<-ggplot(ocdMelt, aes(Outcome_Measure,perc))+
  geom_bar(aes(fill= Frequency),stat="identity")+
  facet_grid(Group~.)+
  scale_y_continuous(labels = scales::percent_format(accuracy = 5L),limits=c(0,1))+
  scale_fill_manual(values=c("grey60","steelblue"),labels=c("Absence", "Presence")) +
  sjPlot::theme_sjplot2() + 
  ylab("Biopsychosocial Compromise")+
  theme(axis.text.x = element_text(size=8,vjust = .5,angle=25), axis.text.y = element_text(size=10),
        axis.title.y = element_text(size = 14)) +
  theme(
    panel.grid.minor = element_blank(), 
    panel.grid.major = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.background = element_blank(),
    line = element_blank())
p5
Figure 5. Bar Plot of Binary Answers By Biopsychosocial Compromise

Figure 5. Bar Plot of Binary Answers By Biopsychosocial Compromise


We generated z scores for those composite scores of multiple variables, in order to make them comparable. In the following figure we may observe the distribution of continuous variables among users with different levels of biopsychosocial compromise.


library(reshape2)  # for melt() function
library(ggplot2)
# First we need to restructure the data into long format:
ocdMelt <- melt(CONS_C1_df_dup_JUL_2020_match_top_sel[c('compromiso_biopsicosocial',"salud_mean_z","dosis_mean_z", "total_mean_z","total_vif_z")], id=c('compromiso_biopsicosocial'))
Warning: attributes are not identical across measure variables; they will be
dropped
names(ocdMelt) <- c('Group', 'Outcome_Measure', 'Frequency')
# plot
ocdBoxplot <- ggplot(ocdMelt, aes(Group, Frequency, color = Outcome_Measure)) + 
  geom_jitter(alpha=0.05) +
  geom_boxplot() + 
  labs(x='Biopsychosocial Compromise', y='Scores', color='Variable') + scale_y_continuous(breaks=seq(0,50, by=5))+  
  sjPlot::theme_sjplot2() + 
  scale_color_manual(values=c("darkslategray4","rosybrown","dimgray","coral3","grey60"),na.value="black") + #"Set2"
  theme(
    panel.grid.minor = element_blank(), 
    panel.grid.major = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.background = element_blank(),
    line = element_blank())
ocdBoxplot
Figure 6. Box Plot of Continuous Answers By Biopsychosocial Compromise

Figure 6. Box Plot of Continuous Answers By Biopsychosocial Compromise


As seen above, users categorized with a severe Biopsychosocial Compromise showed slightly greater patterns of consumption and lower average health in comparison to users with a mild Biopsychosocial Compromise. However, and at least in an univariate level, they seem to not differ much between the different levels.


We checked the multivariate normality of these variables, in order to choose the technique to perform the analyses that may let us infer whether that users in different levels of Biopsychosocial Compromise provene from the same population.


Table 2. Multivariate Normality Test
Test Statistic P Value Normality
Mardia Skewness 49,049.954 0 NO
Mardia Kurtosis 350.767 0 NO
MVN NA NA NO
Note. Using z-scores distributions of Average Health, Average Frequency of Use in the last Month and Average Dose in the las Month


As seen above, variables did not follow a normal distribution.


Non-parametric Multivariate Test

Considering that our variables did not follow a normal distribution, and the test is composed by a mixture of different variable types, we performed a nonparametric test for multivariate samples. The variability of responses on each of the variables of top was compared between the groups using a non-parametric multivariate ANOVA (MANOVA (Burchett et al., 2017)) with biopsychosocial compromise as between-subject factor (discarding inter-observer variability). The procedure generates a randomization test of 5,000 resamples to compute F-statistic (pseudo-F). Additionally, these models have the advantage that allow for differences in between-group variation, is not sensitive to multicollinearity, allows for the inclusion of more variables, and tolerate sparse data.


#mvnormtest::mshapiro.test()

dependent.vars <- cbind(CONS_C1_df_dup_JUL_2020_match_top_sel$salud_mean_z,CONS_C1_df_dup_JUL_2020_match_top_sel$dosis_mean_z, CONS_C1_df_dup_JUL_2020_match_top_sel$total_mean_z)
root.manova <- manova(dependent.vars ~ CONS_C1_df_dup_JUL_2020_match_top_sel$compromiso_biopsicosocial, intercept=F)
summary(root.manova,test="Pillai")
#test='Pillai'). The manova() function provides four multivariate tests by setting the test argument to either Pillai, #Wilks, Hotelling-Lawley, or Roy.

#La prueba de Wilk es la más comúnmente utilizada debido a que fue la primera prueba en ser derivada y a que posee una aproximación F bastante conocida.
#Lawley-Hotelling es también conocida como el estadístico T2 generalizado de Hotelling.
#La prueba de Pillai es la mejor prueba a utilizar en la mayoría de las situaciones. La prueba de Pillai proporciona resultados similares a los de las pruebas de Wilks y de Lawley-Hotelling.
#La prueba de la raíz más grande de Roy es la mejor cuando los vectores medios son colineales. La prueba de Roy no tiene una aproximación F satisfactoria.

summary.aov(root.manova)
CONS_C1_df_dup_JUL_2020_match_top_sel2<-
  CONS_C1_df_dup_JUL_2020_match_top_sel%>%
      dplyr::mutate_at(vars(hurto, robo,  venta_drogas, rina,  otro, lugar_vivir, vivienda,vif), ~as.numeric(as.factor(as.numeric(.))))
#https://stackoverflow.com/questions/59166335/need-help-using-npmv-package-nonparametric-testing-in-r
#https://cran.r-project.org/web/packages/npmv/npmv.pdf
#R.Ellis, A., W.Burchett, W., W.Harrar, S., & Bathke, A. (2017). Nonparametric Inference for Multivariate Data: the Package npmv. JOURNAL OF STATISTICAL SOFTWARE. https://doi.org/10.18637/jss.v076.i04
invisible(
prueba_anormalidad1<-npmv::nonpartest(dosis_mean_z|total_mean_z|salud_mean_z|hurto|robo|venta_drogas|rina|otro|lugar_vivir|vivienda|total_vif_z~compromiso_biopsicosocial, CONS_C1_df_dup_JUL_2020_match_top_sel2,plots=F,permreps=5000,releffects=TRUE)
)
invisible(
prueba_anormalidad2<-npmv::nonpartest(dosis_mean_z|total_mean_z|salud_mean_z|hurto|robo|venta_drogas|rina|otro|lugar_vivir|vivienda|total_vif_z~compromiso_biopsicosocial+origen_ingreso_mod, CONS_C1_df_dup_JUL_2020_match_top_sel2,plots=F,permreps=5000,releffects=TRUE)
)
invisible(
prueba_anormalidad3<-npmv::nonpartest(dosis_mean_z|total_mean_z|salud_mean_z|hurto|robo|venta_drogas|rina|otro|lugar_vivir|vivienda|total_vif_z~compromiso_biopsicosocial+origen_ingreso_mod+dg_trs_fisico, CONS_C1_df_dup_JUL_2020_match_top_sel2,plots=F,permreps=5000,releffects=TRUE)
)
invisible(
prueba_anormalidad4<-npmv::nonpartest(dosis_mean_z|total_mean_z|salud_mean_z|hurto|robo|venta_drogas|rina|otro|lugar_vivir|vivienda|total_vif_z~ compromiso_biopsicosocial+origen_ingreso_mod+dg_trs_fisico+dg_cie_10, CONS_C1_df_dup_JUL_2020_match_top_sel2,plots=F,permreps=5000,releffects=TRUE)
)

#
prueba_anormalidad_res1<-data.table::data.table(prueba_anormalidad1$results, keep.rownames=T)
prueba_anormalidad_res2<-data.table::data.table(prueba_anormalidad2$results, keep.rownames=T)
prueba_anormalidad_res3<-data.table::data.table(prueba_anormalidad3$results, keep.rownames=T)
prueba_anormalidad_res4<-data.table::data.table(prueba_anormalidad4$results, keep.rownames=T)
#INOCRPORË LAS VAR DE CONTROL. SON TODAS DE TIPO FACTOR

sum_df_npmv<-data.frame(cbind(Model=
                                paste0("npmv",c("Raw","Raw","Raw","Raw","2","2","2","2","3","3","3","3","4","4","4","4")),
  rbind(prueba_anormalidad_res1,
      prueba_anormalidad_res2,
      prueba_anormalidad_res3,
      prueba_anormalidad_res4)))%>%
  dplyr::rename("Parameter"="rn")

sum_df_npmv%>% 
  dplyr::mutate(Test.Statistic=round(as.numeric(as.character(Test.Statistic)),3))%>% 
  knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
               caption="Table 3. Analysis of Variance",
                 align =rep('c', 101),
               col.names=c("Model", "Parameter", "Statistic", "DF1", "DF2","p-value", "Permutation Test p-value"))%>%
   kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 8) %>%
  kableExtra::add_footnote( c("Note. Model with average dose (z score), average frequency of use (z score), theft, burglary, drug selling, fights, other, place to live, housing condition and domestic violence. Raw= No control variables; Controlling for Motive of Admission(1); Controlling for Motive of Admission & Physical Condition(2); Controlling for Motive of Admission, Physical Condition & Psychiatric Condition(3)"), notation = "none") %>%
  kableExtra::scroll_box(width = "100%", height = "375px")
Table 3. Analysis of Variance
Model Parameter Statistic DF1 DF2 p-value Permutation Test p-value
npmvRaw ANOVA type test p-value 165.153 10.210 22,736.373 0 0
npmvRaw McKeon approx. for the Lawley Hotelling Test 48.867 22.000 10,909.062 0 0
npmvRaw Muller approx. for the Bartlett-Nanda-Pillai Test 45.850 22.010 12,737.991 0 0
npmvRaw Wilks Lambda 47.360 22.000 12,732.000 0 0
npmv2 ANOVA type test p-value 18.621 23.928 18,350.079 0 0
npmv2 McKeon approx. for the Lawley Hotelling Test 11.546 48.000 18,701.236 0 0
npmv2 Muller approx. for the Bartlett-Nanda-Pillai Test 11.458 48.041 25,481.970 0 0
npmv2 Wilks Lambda 11.507 48.000 24,512.980 0 0
npmv3 ANOVA type test p-value 16.327 13.184 26,222.112 0 0
npmv3 McKeon approx. for the Lawley Hotelling Test 7.570 26.000 11,132.688 0 0
npmv3 Muller approx. for the Bartlett-Nanda-Pillai Test 7.508 26.012 12,733.989 0 0
npmv3 Wilks Lambda 7.539 26.000 12,728.000 0 0
npmv4 ANOVA type test p-value 80.904 7.610 46,227.035 0 0
npmv4 McKeon approx. for the Lawley Hotelling Test 34.028 14.000 6,364.000 0 0
npmv4 Muller approx. for the Bartlett-Nanda-Pillai Test 34.023 14.002 6,363.998 0 0
npmv4 Wilks Lambda 34.028 14.000 6,364.000 0 0
Note. Model with average dose (z score), average frequency of use (z score), theft, burglary, drug selling, fights, other, place to live, housing condition and domestic violence. Raw= No control variables; Controlling for Motive of Admission(1); Controlling for Motive of Admission & Physical Condition(2); Controlling for Motive of Admission, Physical Condition & Psychiatric Condition(3)
##NO OCUPAR: MUY COMPLICADO
#The function ssnonpartest performs an all-subset algorithm to determine which variables cause significant effects, and between which factor levels. See the examples below for some basic uses and look in the help pages for each function for a much more detailed look.
#invisible(
#npmv::ssnonpartest(dosis_mean|total_mean|trabajo_edu_mean~compromiso_biopsicosocial,CONS_C1_df_dup_JUL_2020_match_top_sel,test=c(1,0,0,0),alpha=.05,factors.and.variables=TRUE)
#)



#https://www.rdocumentation.org/packages/vegan/versions/2.4-2/topics/adonis
#t looks like the jaccard distance is really only useful for binary data (presence/absence) while The bray-curtis matrix has been found to be robust for many abundance type data sets, especially those with many paired zeros like stomach content data can have.

vegan_1<- 
  vegan::adonis(CONS_C1_df_dup_JUL_2020_match_top_sel2[,c("dosis_mean_z", "total_mean_z", "salud_mean_z", "hurto", "robo", "venta_drogas", "rina", "otro", "lugar_vivir", "vivienda", "total_vif_z"),drop = FALSE] ~ compromiso_biopsicosocial, data=CONS_C1_df_dup_JUL_2020_match_top_sel,permutations=1000, parallel = 4)
Warning in vegdist(lhs, method = method, ...): results may be meaningless
because data have negative entries in method "bray"
vegan_2<- 
  vegan::adonis(CONS_C1_df_dup_JUL_2020_match_top_sel2[,c("dosis_mean_z", "total_mean_z", "salud_mean_z", "hurto", "robo", "venta_drogas", "rina", "otro", "lugar_vivir", "vivienda", "total_vif_z"),drop = FALSE] ~ compromiso_biopsicosocial+origen_ingreso_mod, data=CONS_C1_df_dup_JUL_2020_match_top_sel,permutations=1000, parallel = 4)
Warning in vegdist(lhs, method = method, ...): results may be meaningless
because data have negative entries in method "bray"
vegan_3<- 
  vegan::adonis(CONS_C1_df_dup_JUL_2020_match_top_sel2[,c("dosis_mean_z", "total_mean_z", "salud_mean_z", "hurto", "robo", "venta_drogas", "rina", "otro", "lugar_vivir", "vivienda", "total_vif_z"),drop = FALSE] ~ compromiso_biopsicosocial+origen_ingreso_mod+dg_trs_fisico, data=CONS_C1_df_dup_JUL_2020_match_top_sel, permutations=1000, parallel = 4)
Warning in vegdist(lhs, method = method, ...): results may be meaningless
because data have negative entries in method "bray"
vegan_4<- 
  vegan::adonis(CONS_C1_df_dup_JUL_2020_match_top_sel2[,c("dosis_mean_z", "total_mean_z", "salud_mean_z", "hurto", "robo", "venta_drogas", "rina", "otro", "lugar_vivir", "vivienda", "total_vif_z"), drop = FALSE] ~ compromiso_biopsicosocial+origen_ingreso_mod+dg_trs_fisico+dg_cie_10,  data=CONS_C1_df_dup_JUL_2020_match_top_sel, permutations=1000, parallel = 4)
Warning in vegdist(lhs, method = method, ...): results may be meaningless
because data have negative entries in method "bray"
sum_df_vegan<-data.frame(cbind(Model=paste0("vegan",c("Raw","Raw","Raw",
                                                      "2","2","2","2",
                                                      "3","3","3","3","3",
                                                      "4","4","4","4","4","4")),
          rbind(data.table::data.table(vegan_1$ aov.tab, keep.rownames=T),
                data.table::data.table(vegan_2$ aov.tab, keep.rownames=T),
                data.table::data.table(vegan_3$ aov.tab, keep.rownames=T),
                data.table::data.table(vegan_4$ aov.tab, keep.rownames=T))))%>%
            dplyr::rename("Parameter"="rn")

sum_df_vegan%>% 
  dplyr::mutate(R2=round(as.numeric(as.character(R2)),3))%>% 
  knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
               caption="Table 4. Analysis of Variance using bray-curtis matrix",
                 align =rep('c', 101),
               col.names=c("Model","Parameter", "Df","SS", "Mean Squares","Pseudo F","R2", "p-value"))%>%
   kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 8) %>%
  kableExtra::add_footnote( c("Note. Model with average dose (z score), average frequency of use (z score), theft, burglary, drug selling, fights, other, place to live, housing condition and domestic violence. Raw= No control variables; Controlling for Motive of Admission(1); Controlling for Motive of Admission & Physical Condition(2); Controlling for Motive of Admission, Physical Condition & Psychiatric Condition(3)"), notation = "none") %>%
  kableExtra::scroll_box(width = "100%", height = "375px")
Table 4. Analysis of Variance using bray-curtis matrix
Model Parameter Df SS Mean Squares Pseudo F R2 p-value
veganRaw compromiso_biopsicosocial 2 11.4921579 5.7460790 172.839357 0.051 0.000999
veganRaw Residuals 6,376 211.9713947 0.0332452 NA 0.949 NA
veganRaw Total 6,378 223.4635527 NA NA 1.000 NA
vegan2 compromiso_biopsicosocial 2 11.4921579 5.7460790 174.720466 0.051 0.000999
vegan2 origen_ingreso_mod 4 2.4137164 0.6034291 18.348410 0.011 0.000999
vegan2 Residuals 6,372 209.5576783 0.0328873 NA 0.938 NA
vegan2 Total 6,378 223.4635527 NA NA 1.000 NA
vegan3 compromiso_biopsicosocial 2 11.4921579 5.7460790 175.033627 0.051 0.000999
vegan3 origen_ingreso_mod 4 2.4137164 0.6034291 18.381297 0.011 0.000999
vegan3 dg_trs_fisico 2 0.4405869 0.2202934 6.710447 0.002 0.000999
vegan3 Residuals 6,370 209.1170914 0.0328284 NA 0.936 NA
vegan3 Total 6,378 223.4635527 NA NA 1.000 NA
vegan4 compromiso_biopsicosocial 2 11.4921579 5.7460790 175.331041 0.051 0.000999
vegan4 origen_ingreso_mod 4 2.4137164 0.6034291 18.412530 0.011 0.000999
vegan4 dg_trs_fisico 2 0.4405869 0.2202934 6.721849 0.002 0.000999
vegan4 dg_cie_10 1 0.3874970 0.3874970 11.823759 0.002 0.000999
vegan4 Residuals 6,369 208.7295944 0.0327727 NA 0.934 NA
vegan4 Total 6,378 223.4635527 NA NA 1.000 NA
Note. Model with average dose (z score), average frequency of use (z score), theft, burglary, drug selling, fights, other, place to live, housing condition and domestic violence. Raw= No control variables; Controlling for Motive of Admission(1); Controlling for Motive of Admission & Physical Condition(2); Controlling for Motive of Admission, Physical Condition & Psychiatric Condition(3)
#dg_cons_dep", "sin_otros_prob_sm","abandono_temp",,"dg_cie_10"
#https://stackoverflow.com/questions/16638297/vegdist-error-message-in-vegan-package

Relationship with outcomes

library(glmulti)
#sin_otros_prob_sm+ dg_cons_dep
res <- glmulti(abandono_temp ~ compromiso_biopsicosocial + origen_ingreso_mod + dg_cie_10 + dg_trs_fisico, #tal vez  no debería meter esta última variable.
               data=CONS_C1_df_dup_JUL_2020_match_top_sel2%>%dplyr::mutate(compromiso_biopsicosocial=case_when(compromiso_biopsicosocial=="2-Moderado"~"Moderate",compromiso_biopsicosocial=="1-Leve"~"Mild",TRUE~"Severe")),
               level=2, #The level of interaction between explanatory variables to be considered. Currently only 1 (only main effects) or 2 (main effects plus all pairwise interactions) are supported.
               fitfunction=glm, crit="aicc", 
               family=binomial,
               marginality=T,#Whether to use the general marginality rule or not. Default is FALSE. With TRUE, all predictors found in interaction terms are also included as main effects.
               confsetsize = 1000 #How many models should be returned in the confidence set of models?
               )
Initialization...
TASK: Exhaustive screening of candidate set.
Fitting...

After 50 models:
Best model: abandono_temp~1+origen_ingreso_mod+dg_cie_10+dg_trs_fisico
Crit= 5789.70750899756
Mean crit= 5942.2170841595

After 100 models:
Best model: abandono_temp~1+origen_ingreso_mod+dg_cie_10+dg_trs_fisico+compromiso_biopsicosocial
Crit= 5779.5616968649
Mean crit= 5866.2943217028
Completed.
#t is not iterative (all models are compared), and is fully automatic (no intervention of the user is required). It is primarily intended to be used for exploratory analyses with many candidate explanatory variables (say 10 to 20)
plot(res)
##########################################################################################
##################DEBIESE INCORPORAR EL TIPO DE PLAN!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
##########################################################################################

#http://www.metafor-project.org/doku.php/tips:model_selection_with_glmulti_and_mumin

#But let's take a look at the top 10 models:
top <- weightable(res) #Prepares a table with model formulas, IC values and IC relative supports
top <- top[top$aicc <= min(top$aicc) + 5,]


############################################################################################

###########################################################################################
plot.glmulti<-function(x, type="p", ...) 
{
  if ( type == "p") {
    plot(x@crits,xlab="Best models",
         ylab=paste("Support (",x@params$crit,")"),pch=19, main="IC profile", ...)
    abline(h=x@crits[1]+2,col="red")
        return(assign("ic_glmulti",data.table::data.table(data.frame(imp),keep.rownames = T),envir = .GlobalEnv))

  } else if (type == "w") {
    ww = exp(-(x@crits - x@crits[1])/2)
    ww=ww/sum(ww)
    plot(ww,xlab="Best models",
         ylab=paste("Evidence weight (",x@params$crit,")"),pch=19, main="Profile of model weights", ...)
    cucu=function(i) sum(ww[1:i])
    wwc=lapply(1:length(ww),cucu)
    abline(v=min(which(wwc>=0.95)),col="red")
  } else if (type=="r") {
    if (length(x@objects)) {
      # shows some diagnostics of the fit
      dev.new()
      par(mfrow=c(2,min(length(x@crits), 5)))
      for (k in 1:min(length(x@crits), 5)) 
        plot(x@objects[[k]],which=c(1), main=deparse(x@formulas[[k]]),...)
      for (k in 1:min(length(x@crits), 5)) 
        plot(x@objects[[k]],which=c(2),...)
    } else warning("Unavailable: use includeobjects=T when calling glmulti.")       
  } else if (type=="s") {
    # plots variable (i.e. terms) importances
    ww = exp(-(x@crits - x@crits[1])/2)
    ww=ww/sum(ww)
    # handle synonymies for interactions
    
    # this translates to unique notations (e.g. x:y and y:x are the same)
    clartou=function(x) {
      sort(strsplit(x, ":")[[1]])-> pieces
      if (length(pieces)>1) paste(pieces[1],":",pieces[2],sep="")
      else x
    }
    # list terms in models
    tet = lapply(x@formulas, function(x) sapply(attr(delete.response(terms(x)),"term.labels"), clartou))
    # all unique terms
    unique(unlist(tet))-> allt
    # importances
    sapply(allt, function(x) sum(ww[sapply(tet, function(t) x%in%t)]))-> imp
    # draw
    #par(oma=c(0,13,0,0))
    barplot(sort(imp),xlab="",xlim=c(0,1), ylab="",horiz=T,las=2, names.arg=allt[order(imp)],main="Model-averaged importance of terms", ...)
    abline(v=0.8, col="red")
    return(assign("var_imp_glmulti",data.table::data.table(data.frame(imp),keep.rownames = T),envir = .GlobalEnv))
  } else    warning("plot: Invalid type argument for plotting glmulti objects.")
}
plot(res, type="s")
############################################################################################
data.frame(res@crits)%>% dplyr::mutate(rn=row_number())%>%
    ggplot(aes(x=rn,y=res.crits))+
    geom_point(color="steelblue")+
   sjPlot::theme_sjplot2() + 
  ylab("Support (aicc)")+
  theme(axis.text.x = element_text(size=8,vjust = .5,angle=25), axis.text.y = element_text(size=10),
        axis.title.y = element_text(size = 14)) +
  theme(
    panel.grid.minor = element_blank(), 
    panel.grid.major = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.background = element_blank(),
    line = element_blank())+
  geom_hline(yintercept=res@crits[1]+2, color="red", linetype=2,size=1)+
  xlab("Best models")
Figure 7. Ordered Models and AICC

Figure 7. Ordered Models and AICC


We selected the model with lower AIC. Of them, the first selected did not contain biopsychosocial compromise (AICc= 5779.53, IC Rel. Support= 0.144). The second contained less interaction terms (did not contain the interaction between cause of admission and biopsychosocial compromise) (AICc= 5780.21, IC Rel. Support= 0.101).


#par(oma=c(0,0,0,24))
ggplot(var_imp_glmulti, aes(x=reorder(rn, -imp),y=imp))+
  geom_bar(stat="identity",alpha=.7, fill="steelblue")+
    sjPlot::theme_sjplot2() + 
  ylab("Model-averaged importance")+
  theme(axis.text.x = element_text(size=6.5,vjust = .5,angle=25), 
        axis.text.y = element_text(size=10),
        axis.title.y = element_text(size = 14)) +
  theme(
    panel.grid.minor = element_blank(), 
    panel.grid.major = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.background = element_blank(),
    line = element_blank())+
  #geom_hline(yintercept=.5, color="red", linetype=2)+
  xlab("Terms")
Figure 8. Model-averaged variable importance of terms

Figure 8. Model-averaged variable importance of terms

df_coef_glmulti<- data.frame()
for (i in 1:nrow(summary(res@objects[[2]])$coefficients)){
    df<-cbind(attr(summary(res@objects[[2]])$coefficients,"dimnames")[[1]][i],
            t(data.frame(exp(summary(res@objects[[2]])$coefficients[i,1] + 
          +     qnorm(c(0.025,.5,0.975)) * summary(res@objects[[2]])$coefficients[i,2]))))
    df_coef_glmulti<-rbind(df_coef_glmulti,df)
          }
df_coef_glmulti<-data.table::data.table(df_coef_glmulti)

data.table::data.table(summary(res@objects[[2]])$coefficients,keep.rownames = T)%>%
  dplyr::left_join(df_coef_glmulti,by=c("rn"="V1"))%>%
  dplyr::select(-V3)%>%
  dplyr::mutate(V2=round(as.numeric(as.character(V2)),2))%>%
  dplyr::mutate(V4=round(as.numeric(as.character(V4)),2))%>%
  dplyr::mutate(Estimate=round(as.numeric(as.character(Estimate)),2))%>%
  dplyr::mutate(OR=round(exp(as.numeric(as.character(Estimate))),2))%>%
  #exp(summary(m)$coefficients["DSH",1] + 
#+     qnorm(c(0.025,0.5,0.975)) * summary(m)$coefficients["DSH",2])
#[1]  1.472098  4.197368 11.967884
  dplyr::mutate(`Std. Error`=round(as.numeric(as.character(`Std. Error`)),2))%>% 
  dplyr::mutate(`z value`=round(as.numeric(as.character(`z value`)),2))%>% 
  dplyr::mutate(`Pr(>|z|)`=round(as.numeric(as.character(`Pr(>|z|)`)),4))%>% 
  dplyr::rename("Terms"="rn", "OR_lo"="V2","OR_up"="V4")%>%
  dplyr::select(Terms,Estimate,OR,OR_lo,OR_up, everything())%>%
  knitr::kable(.,format = "html", format.args = list(decimal.mark = ".", big.mark = ","),
               caption="Table 6. Logistic Regression Coefficients",
                 align =rep('c', 101))%>%
   kableExtra::kable_styling(bootstrap_options = c("striped", "hover"),font_size = 8)%>%
  kableExtra::add_footnote( c(paste("Note. AICc=",round(res@objects[[2]]$aic,3)),
                              paste("Defined model=",top$model[2])), notation = "none")%>%
  kableExtra::scroll_box(width = "100%", height = "375px")
Table 6. Logistic Regression Coefficients
Terms Estimate OR OR_lo OR_up Std. Error z value Pr(>|z|)
(Intercept) -2.17 0.11 0.08 0.16 0.17 -13.03 0.0000
origen_ingreso_modDerivación Asistida -0.23 0.79 0.62 1.01 0.12 -1.92 0.0554
origen_ingreso_modOtros -0.50 0.61 0.43 0.86 0.18 -2.77 0.0057
origen_ingreso_modSector Justicia -0.34 0.71 0.55 0.92 0.13 -2.61 0.0090
origen_ingreso_modSector Salud -0.36 0.70 0.60 0.82 0.08 -4.37 0.0000
dg_cie_101 1.49 4.44 3.76 5.27 0.09 17.34 0.0000
dg_trs_fisico1 -0.68 0.51 0.33 0.79 0.22 -3.04 0.0023
dg_trs_fisico2 -0.81 0.44 0.18 1.12 0.47 -1.72 0.0846
compromiso_biopsicosocialModerate 0.03 1.03 0.73 1.45 0.17 0.19 0.8520
compromiso_biopsicosocialSevere 0.08 1.08 0.75 1.59 0.19 0.44 0.6626
dg_trs_fisico1:compromiso_biopsicosocialModerate 0.29 1.34 0.84 2.16 0.24 1.22 0.2242
dg_trs_fisico2:compromiso_biopsicosocialModerate 0.24 1.27 0.47 3.37 0.50 0.47 0.6378
dg_trs_fisico1:compromiso_biopsicosocialSevere 0.57 1.77 1.07 2.93 0.26 2.23 0.0258
dg_trs_fisico2:compromiso_biopsicosocialSevere 0.13 1.14 0.42 3.08 0.51 0.26 0.7920
Note. AICc= 5780.205
Defined model= abandono_temp ~ 1 + origen_ingreso_mod + dg_cie_10 + dg_trs_fisico + compromiso_biopsicosocial + dg_trs_fisico:compromiso_biopsicosocial
#abandono_temp ~ 1 + compromiso_biopsicosocial + origen_ingreso_mod + dg_cie_10 + dg_trs_fisico + dg_trs_fisico:compromiso_biopsicosocial
#print(res@formulas[[2]])


The next figures show the marginal effects of the interaction terms among the variables.


library(lsmeans)
refR_1 <- lsmeans(res@objects[[2]], specs = c("compromiso_biopsicosocial","dg_trs_fisico"))
ref_dfR_1 <- as.data.frame(cbind(summary(refR_1)[1:2],exp(summary(refR_1)[,c(3,6,7)])))

g4R_1 <- ggplot(ref_dfR_1, aes(x=compromiso_biopsicosocial, y=lsmean,group=dg_trs_fisico, colour=dg_trs_fisico))+
geom_errorbar(aes(ymin=asymp.LCL, ymax=asymp.UCL), width=.1,position=position_dodge(0.1), size=1) +
geom_line(position=position_dodge(0.1), size=.5)+
geom_point(position=position_dodge(0.1), size=2)+
  xlab("Biopsychosocial Compromise")+
  ylab("Odds Ratio (OR) of Early Drop-out")+
  sjPlot::theme_sjplot2() +
  geom_rect_interactive(alpha = 0.1, xmin=.1, xmax=.1, ymin=.1,ymax=.1) +
    # Remove plot elements added by geom_rect_interactive
  theme(legend.position="bottom")+
  guides(color=guide_legend(ncol=4))+
  labs(color="Cause of Discharge")+
  scale_colour_brewer(name="Cumulative Percentage", palette = "Set1",labels=c("No Dg.", "In study", "Dg."))+
  theme(legend.title = element_blank())
  #theme(axis.text = element_blank(), axis.ticks = element_blank())

ggiraph(code = print(g4R_1))

Figure 9. Interactive Effects in the Prob. of Early Drop-outv with Diagnose of Physical Diagnosis

library(lsmeans)
refR_2 <- lsmeans(res@objects[[2]], specs = c("dg_cie_10","origen_ingreso_mod"))
ref_dfR_2 <- as.data.frame(cbind(summary(refR_2)[1:2],exp(summary(refR_2)[,c(3,6,7)])))

g4R_2 <- ggplot(ref_dfR_2, aes(x=dg_cie_10, y=lsmean,group=origen_ingreso_mod, colour=origen_ingreso_mod))+
geom_errorbar(aes(ymin=asymp.LCL, ymax=asymp.UCL), width=.1,position=position_dodge(0.1), size=1) +
geom_line(position=position_dodge(0.1), size=.5)+
geom_point(position=position_dodge(0.1), size=2)+
  xlab("Psychiatric Condition")+
  ylab("Odds Ratio (OR) of Early Drop-out")+
  sjPlot::theme_sjplot2() +
  geom_rect_interactive(alpha = 0.1, xmin=.1, xmax=.1, ymin=.1,ymax=.1) +
    # Remove plot elements added by geom_rect_interactive
  theme(legend.position="bottom")+
  guides(color=guide_legend(ncol=4,name = "Cause of Discharge"))+
  labs(color="Motive of Admission to Treatment")+
  scale_colour_brewer(palette = "Set1")+#,labels=c("No Dg.", "In study", "Dg."))+
  theme(legend.title = element_blank())
  #theme(axis.text = element_blank(), axis.ticks = element_blank())

ggiraph(code = print(g4R_2))
library(lsmeans)
refR_3 <- lsmeans(res@objects[[2]], specs = c("dg_cie_10","dg_trs_fisico"))
ref_dfR_3 <- as.data.frame(cbind(summary(refR_3)[1:2],exp(summary(refR_3)[,c(3,6,7)])))

g4R_3 <- ggplot(ref_dfR_3, aes(x=dg_cie_10, y=lsmean,group=dg_trs_fisico, colour=dg_trs_fisico))+
geom_errorbar(aes(ymin=asymp.LCL, ymax=asymp.UCL), width=.1,position=position_dodge(0.1), size=1) +
geom_line(position=position_dodge(0.1), size=.5)+
geom_point(position=position_dodge(0.1), size=2)+
  xlab("Physical Condition")+
  ylab("Odds Ratio (OR) of Early Drop-out")+
  sjPlot::theme_sjplot2() +
  geom_rect_interactive(alpha = 0.1, xmin=.1, xmax=.1, ymin=.1,ymax=.1) +
    # Remove plot elements added by geom_rect_interactive
  theme(legend.position="bottom")+
  guides(color=guide_legend(ncol=4,name = "Cause of Discharge"))+
  labs(color="Motive of Admission to Treatment")+
  scale_colour_brewer(palette = "Set1",labels=c("No Dg.", "In study", "Dg."))+
  theme(legend.title = element_blank())
  #theme(axis.text = element_blank(), axis.ticks = element_blank())

ggiraph(code = print(g4R_3))
save.image("G:/Mi unidad/Alvacast/SISTRAT 2019 (github)/biopsychosoc_comp.RData")
sessionInfo()
R version 4.0.2 (2020-06-22)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18363)

Matrix products: default

locale:
[1] LC_COLLATE=Spanish_Chile.1252  LC_CTYPE=Spanish_Chile.1252   
[3] LC_MONETARY=Spanish_Chile.1252 LC_NUMERIC=C                  
[5] LC_TIME=Spanish_Chile.1252    

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] gdtools_0.2.2           glmulti_1.0.8           leaps_3.1              
 [4] rJava_0.9-12            reshape2_1.4.4          polycor_0.7-10         
 [7] DiagrammeR_1.0.6.1.9000 MVN_5.8                 citr_0.3.2             
[10] ggcorrplot_0.1.3        psych_2.0.7             ggiraph_0.7.0          
[13] finalfit_1.0.1          lsmeans_2.30-0          emmeans_1.4.7          
[16] choroplethrAdmin1_1.1.1 choroplethrMaps_1.0.1   choroplethr_3.6.3      
[19] acs_2.1.4               XML_3.99-0.3            RColorBrewer_1.1-2     
[22] panelr_0.7.3            lme4_1.1-23             Matrix_1.2-18          
[25] dplyr_1.0.0             data.table_1.12.8       codebook_0.9.2         
[28] Statamarkdown_0.4.5     devtools_2.3.0          usethis_1.6.1          
[31] sqldf_0.4-11            RSQLite_2.2.0           gsubfn_0.7             
[34] proto_1.0.0             broom_0.7.0             zoo_1.8-8              
[37] altair_4.0.1            rbokeh_0.5.1            janitor_2.0.1          
[40] plotly_4.9.2.1          kableExtra_1.1.0        Hmisc_4.4-0            
[43] Formula_1.2-3           survival_3.1-12         lattice_0.20-41        
[46] ggplot2_3.3.1           stringr_1.4.0           stringi_1.4.6          
[49] tidyr_1.1.0             knitr_1.29              matrixStats_0.56.0     
[52] boot_1.3-25            

loaded via a namespace (and not attached):
  [1] vcd_1.4-7             class_7.3-17          ps_1.3.3             
  [4] lmtest_0.9-37         rprojroot_1.3-2       crayon_1.3.4         
  [7] laeken_0.5.1          MASS_7.3-51.6         nlme_3.1-148         
 [10] backports_1.1.8       rlang_0.4.6           readxl_1.3.1         
 [13] performance_0.4.6     nloptr_1.2.2.1        callr_3.4.3          
 [16] rjson_0.2.20          ggmap_3.0.0           bit64_0.9-7          
 [19] glue_1.4.1            sjPlot_2.8.4          processx_3.4.2       
 [22] classInt_0.4-3        tcltk_4.0.2           haven_2.3.1          
 [25] tidyselect_1.1.0      NADA_1.6-1.1          rio_0.5.16           
 [28] sf_0.9-3              sjmisc_2.8.5          chron_2.3-55         
 [31] xtable_1.8-4          magrittr_1.5          evaluate_0.14        
 [34] RgoogleMaps_1.4.5.3   cli_2.0.2             rstudioapi_0.11      
 [37] miniUI_0.1.1.1        sp_1.4-2              robCompositions_2.2.1
 [40] rpart_4.1-15          jtools_2.0.5          pls_2.7-2            
 [43] zCompositions_1.3.4   sjlabelled_1.1.5      RJSONIO_1.3-1.4      
 [46] maps_3.3.0            shiny_1.5.0           gistr_0.5.0          
 [49] xfun_0.16             parameters_0.7.0      pkgbuild_1.1.0       
 [52] cluster_2.1.0         sgeostat_1.0-27       tibble_3.0.1         
 [55] permute_0.9-5         png_0.1-7             reshape_0.8.8        
 [58] withr_2.2.0           bitops_1.0-6          ranger_0.12.1        
 [61] plyr_1.8.6            cellranger_1.1.0      pcaPP_1.9-73         
 [64] e1071_1.7-3           coda_0.19-3           pillar_1.4.6         
 [67] mvoutlier_2.0.9       multcomp_1.4-13       fs_1.4.1             
 [70] flexmix_2.3-15        kernlab_0.9-29        vctrs_0.3.1          
 [73] ellipsis_0.3.1        generics_0.0.2        nortest_1.0-4        
 [76] rgdal_1.5-8           tools_4.0.2           foreign_0.8-80       
 [79] munsell_0.5.0         fastmap_1.0.1         compiler_4.0.2       
 [82] pkgload_1.1.0         abind_1.4-5           httpuv_1.5.4         
 [85] tigris_0.9.4          npmv_2.4.0            sessioninfo_1.1.1    
 [88] gridExtra_2.3         visNetwork_2.0.9      later_1.1.0.1        
 [91] jsonlite_1.6.1        WDI_2.6.0             GGally_1.5.0         
 [94] scales_1.1.1          carData_3.0-4         estimability_1.3     
 [97] lazyeval_0.2.2        promises_1.1.1        car_3.0-8            
[100] latticeExtra_0.6-29   reticulate_1.16       effectsize_0.3.1     
[103] checkmate_2.0.0       rmarkdown_2.3         openxlsx_4.1.5       
[106] sandwich_2.5-1        moments_0.14          statmod_1.4.34       
[109] webshot_0.5.2         forcats_0.5.0         pander_0.6.3         
[112] yaml_2.2.1            systemfonts_0.2.3     prabclus_2.3-2       
[115] htmltools_0.5.0       memoise_1.1.0         modeltools_0.2-23    
[118] viridisLite_0.3.0     digest_0.6.25         rrcov_1.5-2          
[121] assertthat_0.2.1      mime_0.9              rappdirs_0.3.1       
[124] repr_1.1.0            bayestestR_0.6.0      units_0.6-6          
[127] remotes_2.1.1         vegan_2.5-6           blob_1.2.1           
[130] splines_4.0.2         labeling_0.3          fpc_2.2-7            
[133] cvTools_0.3.2         hms_0.5.3             modelr_0.1.8         
[136] colorspace_1.4-1      base64enc_0.1-3       mnormt_2.0.0         
[139] tmvnsim_1.0-2         nnet_7.3-14           Rcpp_1.0.4.6         
[142] mclust_5.4.6          mvtnorm_1.1-1         fansi_0.4.1          
[145] VIM_6.0.0             truncnorm_1.0-8       R6_2.4.1             
[148] grid_4.0.2            lifecycle_0.2.0       acepack_1.4.1        
[151] labelled_2.5.0        zip_2.0.4             curl_4.3             
[154] pryr_0.1.4            minqa_1.2.4           testthat_2.3.2       
[157] snakecase_0.11.0      robustbase_0.93-6     skimr_2.1.1          
[160] desc_1.2.0            TH.data_1.0-10        htmlwidgets_1.5.1    
[163] purrr_0.3.4           crosstalk_1.1.0.1     mgcv_1.8-31          
[166] rvest_0.3.5           insight_0.8.4         htmlTable_1.13.3     
[169] codetools_0.2-16      lubridate_1.7.9       prettyunits_1.1.1    
[172] vegawidget_0.3.1      gtable_0.3.0          DBI_1.1.0            
[175] stats4_4.0.2          httr_1.4.2            highr_0.8            
[178] KernSmooth_2.23-17    farver_2.0.3          uuid_0.1-4           
[181] diptest_0.75-7        hexbin_1.28.1         mice_3.9.0           
[184] xml2_1.3.2            ggeffects_0.14.3      readr_1.3.1          
[187] sROC_0.1-2            energy_1.7-7          DEoptimR_1.0-8       
[190] bit_1.1-15.2          sjstats_0.18.0        jpeg_0.1-8.1         
[193] pkgconfig_2.0.3       maptools_1.0-1       

References

1. Revelle W. psych: Procedures for personality and psychological research.

2. Fox J. polycor: Polychoric and Polyserial Correlations.

3. Silva J da, Ventura CAA, Vargens OM da C, et al. Illicit drug use in seven Latin American countries: critical perspectives of families and familiars. Revista Latino-Americana de Enfermagem 2009; 17: 763–769.

4. Basu D, Ghosh A, Sarkar S, et al. Initial treatment dropout in patients with substance use disorders attending a tertiary care de-addiction centre in north India. The Indian journal of medical research 2017; 146: S77–S84.

5. Lappan SN, Brown AW, Hendricks PS. Dropout rates of in-person psychosocial substance use disorder treatments: a systematic review and meta-analysis. Addiction 2020; 115: 201–217.

6. Pérez-Gómez A, Mejía-Trujillo J. The Evolution of Alcohol and Drug Prevention Strategies in Latin America. In: Romano JL, Israelashvili M (eds) The cambridge handbook of international prevention science. Cambridge: Cambridge University Press, pp. 753–779.

7. Ronzani TM, Fuentes-Mejía C, Mota DCB, et al. Intervenciones breves para el abuso de sustancias en América Latina: una revisión sistemática. Psicologia em Estudo; 24, http://www.scielo.br/scielo.php?script=sci{\_}arttext{\&}pid=S1413-73722019000100232{\&}nrm=iso (2019).

8. Heinze G, Armas-Castañeda G. Public policies on the use of drugs in Mexico and Latin America. Drug Science, Policy and Law 2015; 2: 2050324515611587.

9. Calcagno V, Mazancourt C de. glmulti: an R package for easy automated model selection with (generalized) linear models. Epub ahead of print 2010. DOI: 10.18637/jss.v034.i12.

10. Dixon P. VEGAN, a package of R functions for community ecology. Journal of Vegetation Science 2003; 14: 927–930.

11. Team RC. R: A language and environment for statistical computing., http://www.r-project.org/ (2020).

12. Inter-American Drug Abuse Control Commission [CICAD]. Report on Drug Use in the Americas 2019. 2019. OEA/Ser.L/XIV.6.6, Washington, D.C.: Organization of American States (OAS), http://www.cicad.oas.org/main/pubs/Report on Drug Use in the Americas 2019.pdf (2019).

13. Korkmaz D. S, Zararsiz G. MVN: An R Package for Assessing Multivariate Normality. The R Journal 2014; 6: 151–162.

14. Burchett W, Ellis A, Harrar S, et al. Nonparametric Inference for Multivariate Data: The R Package npmv. Journal of Statistical Software; 76. Epub ahead of print 2017. DOI: 10.18637/jss.v076.i04.

15. Soto-Brandt G, Portilla Huidobro R, Huepe Artigas D, et al. [Validity evidence of the Alcohol, Smoking and Substance Involvement Screening Test (ASSIST) in Chile]. Adicciones 2014; 26: 291–302.

16. Hansen Z. A Phenomenological Investigation of Clinical Intuition among Alcohol and Drug Counselors. 2015.

17. Stone J, Marsh A, Dale A, et al. Counselling Guideleines: Fourth Edition: Alchol and Other Drug Issues. Western Australian Government - Mental Health Commission, https://books.google.cl/books?id=K{\_}k2xQEACAAJ (2019).

18. Group CG on DM, Working DU2IE. Drug misuse and dependence: UK guidelines on clinical management. Public Health England. Alcohol, Drugs & Tobacco Division, https://assets.publishing.service.gov.uk/government/uploads/system/uploads/attachment{\_}data/file/673978/clinical{\_}guidelines{\_}2017.pdf (2017).

19. Center for Substance Abuse T. SAMHSA/CSAT Treatment Improvement Protocols. In: Substance abuse treatment for persons with co-occurring disorders. Rockville (MD): Substance Abuse; Mental Health Services Administration (US), 2005.

20. Maunder RG, Wiesenfeld L, Lawson A, et al. The Relationship Between Childhood Adversity and Other Aspects of Clinical Complexity in Psychiatric Outpatients. Journal of Interpersonal Violence 2019; 0886260519865968.

21. Hsieh M-H, Tsai S-L, Tsai C-H, et al. What Is the Addiction World Like? Understanding the Lived Experience of the Individuals’ Illicit Drug Addiction in Taiwan. Perspectives in Psychiatric Care 2017; 53: 47–54.

22. Messas G, Fukuda L, Pienkos E. A Phenomenological Contribution to Substance Misuse Treatment: Principles for Person-Centered Care. Psychopathology 2019; 52: 1–9.

23. Copoeru I. Understanding Addiction: A Threefold Phenomenological Approach. Human Studies 2014; 37: 335–349.

24. Pachado MP, Scherer JN, Guimarães LSP, et al. Markers for Severity of Problems in Interpersonal Relationships of Crack Cocaine Users from a Brazilian Multicenter Study. Psychiatric Quarterly 2018; 89: 923–936.

25. Ministerio de Salud[MINSAL]. Norma y orientaciones técnicas de los planes de tratamiento y rehabilitación para personas adultas con problemas derivados del consumo de drogas. 391, https://www.senda.gob.cl/wp-content/uploads/2012/08/OrientacionesTecnicas{\_}CentrosdeTratamiento.pdf (2012).

26. Ministerio de Salud[MINSAL]. Norma técnica para el tratamiento integral de adolescentes infractores de ley con consumo problemático de alcohol - drogas y otros trastornos de salud mental. 391, https://www.minsal.cl/portal/url/item/71e5abf67b425395e04001011f017d2e.pdf (2006).

27. Valdivieso M. Sistematización de las prácticas de un programa de rehabilitación en drogas y alcohol para mujeres, con enfoque de género, que ha venido realizando el equipo de profesionales y técnicos del CTR del hospital de Peñablanca. PhD thesis, Universidad de Chile, http://repositorio.uchile.cl/handle/2250/135697 (2014).

28. Consejo Nacional para el Control de Estupefacientes[CONACE]. Modelo de intervención en personas con consumo problemático de sustancias psicoactivas, recluidas en los establecimientos penitenciarios chilenos. Tomo II. Programa de Tratamiento, Rehabilitación y Reinserción Social, para internos/as con consumo problemático de sustancias psicoactivas. Ministerio del Interior, http://sistemas.senda.gob.cl/sistema-monitoreo/biblioteca/files/Documentos/ESTRATEGIAS NORMAS ORIENTACIONES/1 Orientaciones y Normas/Nacional/Senda/Modelo de intervenci{\'{o}}n de personas recluidas en establecimientos penitenciarios VOL 2.pdf (2005).

29. Consejo Nacional para el Control de Estupefacientes [CONACE]. Orientaciones Técnicas. Tratamiento del consumo problemático de alcohol y drogas y otros trastornos de salud mental en adolescentes infractores de ley. 2007.

30. Consejo Nacional para el Control de Estupefacientes[CONACE]. Modelo de intervención en personas con consumo problemático de sustancias psicoactivas, recluidas en los establecimientos penitenciarios chilenos. Tomo I. Elementos Teóricos del Programa de Tratamiento, Rehabilitación y Reinserción Social, para Internos/as con consumo problemático de sustancias psicoactivas. Ministerio del Interior, http://sistemas.senda.gob.cl/sistema-monitoreo/biblioteca/files/Documentos/ESTRATEGIAS NORMAS ORIENTACIONES/1 Orientaciones y Normas/Nacional/Senda/Modelo de intervenci{\'{o}}n de personas recluidas en establecimientos penitenciarios VOL 1.pdf (2005).

31. Huidobro R, Lozier M, Oliva M. Informe resultados: Treatment Outcomes Profile (TOP) Chile. 2016.

32. Tiffany ST, Friedman L, Greenfield SF, et al. Beyond drug use: a systematic consideration of other outcomes in evaluations of treatments for substance use disorders. Addiction (Abingdon, England) 2012; 107: 709–718.

33. National Institute on Drug Abuse[NIDA]. Part 1: The Connection Between Substance Use Disorders and Mental Illness. National Institute on Drug Abuse, https://www.drugabuse.gov/publications/research-reports/common-comorbidities-substance-use-disorders/part-1-connection-between-substance-use-disorders-mental-illness.

34. Simoneau H, Brochu S. Addiction severity index profile of persons who reenter treatment for substance use disorders. Substance Abuse 2017; 38: 432–437.

35. Leung J, Peacock A, Colledge S, et al. A Global Meta-analysis of the Prevalence of HIV, Hepatitis C Virus, and Hepatitis B Virus Among People Who Inject Drugs—Do Gender-Based Differences Vary by Country-Level Indicators? The Journal of Infectious Diseases 2019; 220: 78–90.

36. Lozano ÓM, Rojas AJ, Fernández Calderón F. Psychiatric comorbidity and severity of dependence on substance users: how it impacts on their health-related quality of life? Journal of Mental Health 2017; 26: 119–126.

37. Kopak AM, Combs E, Goodman K, et al. Exposure to Violence and Substance Use Treatment Outcomes Among Female Patients. Substance Use & Misuse 2019; 54: 362–372.

38. Medlock MM, Rosmarin DH, Connery HS, et al. Religious coping in patients with severe substance use disorders receiving acute inpatient detoxification. The American journal on addictions 2017; 26: 744–750.

39. Hildebrandt T, Epstein EE, Sysko R, et al. Using Factor Mixture Models to Evaluate the Type A/B Classification of Alcohol Use Disorders in a Heterogeneous Treatment Sample. Alcoholism, clinical and experimental research 2017; 41: 987–997.

40. Akbari H, Roshanpajouh M, Nourijelyani K, et al. Profile of drug users in the residential treatment centers of Tehran, Iran. Health promotion perspectives 2019; 9: 248–254.

41. Saberi Zafarghandi MB, Khanipour H, Ahmadi SM. Typology of Substance Use Disorder Based on Temperament Dimensions, Addiction Severity, and Negative Emotions. Iranian journal of psychiatry 2018; 13: 184–190.

42. Johannessen DA, Nordfjærn T, Geirdal AØ. Change in psychosocial factors connected to coping after inpatient treatment for substance use disorder: a systematic review. Substance abuse treatment, prevention, and policy 2019; 14: 16.

43. Lopes-Rosa R, Kessler FP, Pianca TG, et al. Predictors of early relapse among adolescent crack users. Journal of Addictive Diseases 2017; 36: 136–143.

44. Klein M. Relapse into opiate and crack cocaine misuse: a scoping review. Addiction Research & Theory 2020; 1–19.

45. Schmid M, Michaud L, Bovio N, et al. Prevalence of somatic and psychiatric morbidity across occupations in Switzerland and its correlation with suicide mortality: results from the Swiss National Cohort (1990–2014). BMC Psychiatry 2020; 20: 324.

46. Gmel G, Akre C, Astudillo M, et al. The Swiss Cohort Study on Substance Use Risk Factors - Findings of two Waves. Sucht 2015; 61: 251–262.

47. Hjemsæter AJ, Bramness JG, Drake R, et al. Mortality, cause of death and risk factors in patients with alcohol use disorder alone or poly-substance use disorders: a 19-year prospective cohort study. BMC Psychiatry 2019; 19: 101.

48. Moradinazar M, Najafi F, Jalilian F, et al. Prevalence of drug use, alcohol consumption, cigarette smoking and measure of socioeconomic-related inequalities of drug use among Iranian people: findings from a national survey. Substance abuse treatment, prevention, and policy 2020; 15: 39.

49. Marín-Navarrete M. R-M, Pérez-López A, Horigian VE. Development and evaluation of addiction treatment programs in Latin America. Current opinion in psychiatry 2018; 31: 306–314.

50. Saberi H.Ahmadi, SM. M. Typology of Substance Use Disorder Based on Temperament Dimensions, Addiction Severity, and Negative Emotions. Iran J Psychiatry 2018; 13: 184–190.

51. Portela A, Garcia L, Goldim J. Adolescencia vulnerable: factores biopsicosociales relacionados al uso de drogas. Revista Bioética 2015; 23: 311–319.

52. Papadimitriou G. The "Biopsychosocial Model": 40 years of application in Psychiatry. Psychiatriki 2017; 28: 107–110.

53. Organization[WHO] WH. World Drug Report 2019 (Sales No. E.19.XI.8). Viena, 2020.

54. Teesson M, Degenhardt L, Hall W. Theories of addiction: Causes and maintenance of addiction. In: Addictions. London: Psychology Press, 2002, pp. 33–47.

55. Vos T, Allen C, Arora M, et al. Global, regional, and national incidence, prevalence, and years lived with disability for 310 diseases and injuries, 1990–2015: a systematic analysis for the Global Burden of Disease Study 2015. The Lancet. Epub ahead of print 2016. DOI: 10.1016/S0140-6736(16)31678-6.

56. Zijlmans EAO, Ark LA van der, Tijmstra J, et al. Methods for Estimating Item-Score Reliability. Applied Psychological Measurement. Epub ahead of print 2018. DOI: 10.1177/0146621618758290.

57. Castillo-Carniglia Á, Marín JD, Soto-Brandt G, et al. Adaptation and Validation of the Instrument Treatment Outcomes Profile to the Chilean Population. Journal of Substance Abuse Treatment. Epub ahead of print 2015. DOI: 10.1016/j.jsat.2015.03.002.

58. Collaborators GBD2D, Incidence I, Prevalence. Global, regional, and national incidence, prevalence, and years lived with disability for 354 diseases and injuries for 195 countries and territories, 1990-2017: a systematic analysis for the Global Burden of Disease Study 2017. Lancet (London, England) 2018; 392: 1789–1858.

59. Polychoric T, Correlations P, Fox AJ. Package ‘polycor’. Statistics.

60. Pacurucu-Castillo SF, Ordóñez-Mancheno JM, Hernández-Cruz A, et al. World Opioid and Substance Use Epidemic: A Latin American Perspective. Psychiatric Research and Clinical Practice. Epub ahead of print 2019. DOI: 10.1176/appi.prcp.20180009.

61. De Boni RB, Peratikos MB, Shepherd BE, et al. Is substance use associated with HIV cascade outcomes in Latin America? PLoS ONE. Epub ahead of print 2018. DOI: 10.1371/journal.pone.0194228.